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1 Introduction 



This paper continues an investigation into the use of multiquadric functions to approximate 
scattered data. An initial report on our efforts appears in [FH93], where we reported that 
a phenomenon that seemed particularly interesting was the tendency of the optimization 
process to result in near-repeated knots. The implication is that the method was attempting 
to build a directional derivative of the basis function into the approximation. This paper 
documents our experiences with an algorithm that detects the occurrence of near-repeated 
knots, and when they occur, replaces basis functions corresponding to near-repeated knots 
with a single multiquadric basis function plus appropriate derivatives of the multiquadric 
basis function, at the pertinent knot. 

While it is useful to have the previous paper at hand, for completeness we will briefly 
review the background necessary to make this report somewhat self-contained. We restrict 
our discussion to functions of two independent variables, the methods are easily extendible 
to arbitrary dimensions, and we expect that many of the conclusions will carry over. 

The scattered data approximation problem is easily described and occurs frequently in 
many branches of science. The problem occurs in any discipline where measurements are 
taken at irregularly spaced values of two or more independent variables, and is especially 
prevalent in environmental sciences. We will suppose that triples of data,(xj, y 3 , z 3 ), j = 
1 , • - • , A r are given, assumed to be measurements (perhaps with error) of an underlying 
function z = f(x,y). The function / is to be approximated by a function F(x,y) from the 
given data. A recent survey of such methods is given in [ F N 9 1 ] . 

Multiquadric functions were introduced for interpolation of scattered data by Hardy 
[HA71]; also see [HA90] for a historical survey and many references. The method is one 
of a class of methods known now as "radial basis function methods'* that includes other 
attractive schemes such as thin plate splines [HD7‘2, DU76, and others]. The basic idea 
of such methods is quite simple, and we describe it in some generality; for purposes of 
being definite it is pertinent to note that for the multiquadric method the radial function 
is h(d) = \[{d 2 + r 2 ). In general, suppose a function of one variable, /? ( c/ ) , where d denotes 
distance, is given. 

For interpolation (that is, exact matching of the given data), a basis function, B 3 {x,y) — 
h(dj) is associated with each data point. Here d 3 = yj{{x ~ Xj) 2 + (y — y^) 2 ), the distance 
from (x,y) to (x 3 ,yj). Thus each basis function is a translate of the radial function, h. The 
approximation is a linear combination of the basis functions, along with some polynomial 
terms that may be necessary in some cases, or may be used to assure that the approximation 
method has polynomial precision. Thus, 

N M 

F(x,y) = Y, a jBj(x,y) + '52b J <lj(x,y) ( 1 ) 

j=\ j=i 

where {^} is a set of M polynomials forming a basis for polynomials of degree < m. The co- 
efficients a } and bj are determined by the linear system of equations prescribing interpolation 
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of the data, and exactness for polynomials of degree < m: 



N M 

y v 'i'Vi) “I" > " bj(Jj (x t , Vi ) = Z{ , l — 1 , ‘ , A 

J = 1 J=1 

N 

Q J 9i(*j» yj) = o, i = 1, • • • , M. 

j=i 



( 2 ) 



For multiquadric basis functions, this system of equations is known to have a unique solution 
for distinct (xj ,?/_,) data (see, for example, [MI86]); while m may be taken as zero (no 
polynomial terms), Micchelli’s results show the multiquadric basis function is positive definite 
of order one, and thus a constant term should be included. We have done so in all our work. If 
higher degree polynomial precision is desired, inclusion of those terms imposes no particular 
burden. 

While interpolation theory is important and indicates something about the suitability 
of the class of functions for approximation purposes, our emphasis here is on least squares 
approximation. This implies using fewer basis functions than there are data points. In 
analogy with univariate cubic splines, it is convenient to refer to the points at which the 
radial basis functions are centered as “knots”, as was done in [MF92], and we do so here. 
If a set of knot points, (uk,Vk),k = l, - I\ , with I\ < N have been specified, then the 
problem of fitting a multiquadric function by least squares is similar to that of solving 
the system of equations corresponding to those above in the least squares sense. We give 
the details. Now, let Bk(x,y) denote the radial basis function associated with the point 
(uk, Vk), Bk(x, y) = \J((x — Uk ) 2 + (y — Vk ) 2 + r 2 ). The system of equations, specialized for 
our case, is now of the form 



K 

) ' Ok Bk(x { , y,) + c z t , i 1, • • • , A 

k=\ 

y o,k = o. 

t=i 



( 3 ) 



There is a question of how to treat the last equation, which guarantees precision for con- 
stants. In [MF92] the corresponding constraint equations were imposed exactly, rather than 
approximately, because of physical considerations. While there is not the corresponding 
physical situation here, we have also imposed the last equation as a constraint. This con- 
straint can be used to reduce the size of the system by solving for a k in terms of the other 
dk and substituting into the first set of equations. 

If the knot points are a subset of the data points, then the system of equations (by 
virtue of containing, as a subset, interpolation equations) has full rank, and thus guarantees 
a unique solution of the least squares problem. When the knot locations may differ from the 
data points, the problem of whether the coefficient matrix is of full rank or not is unknown to 
us, although we feel certain that the matrix is of full rank when the knot points are distinct, 
and have encountered no situations that indicate otherwise. 

The impetus behind our original investigation was to obtain surface approximations that 
would be efficient in subsequent applications. That is, we consider it to be acceptable to 
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expend considerable computational resources to obtain the approximation in a preprocessing 
step. Once obtained, the approximation can be evaluated efficiently, so it would be feasible 
to use it numerous times in an application program. 

In the previous paper [FH93] we found that there was a tendency of the algorithm to 
converge on solutions where the knot points were often close to each other, or near-repeated. 
This occurred both in clusters of two and three knots. At that time, no attempt was made to 
made to investigate this behavior further. The outstanding result of that reference was the 
finding that allowing the knots (and r value) to be determined as part of the minimization 
process yielded approximations that were much more efficient than those obtained either by 
determining the knots a priori or adaptively (we outlined a greedy algorithm for sequentially 
determining knot locations that worked well). While we found that using variable r values 
also yielded good improvement, we do not pursue that idea here. 

Because the greedy algorithm we developed previously was generally used as the starting 
point for the work reported here, we give a brief outline of it. 

a) Obtain the least squares fit by a constant function, the average of the data values. The 

two data points having maximum positive and maximum negative error are taken to 
be the first two knots, (ui,i>j) and (u 2 , i’ 2 ). The knot counter K is set to 2. 

b) The least squares multiquadric fit with I\ knots is obtained, and the residuals are 

computed. 

c) The maximum absolute value of the residuals is found and the location of this residual, 

subject to the minimum knot separation value, is taken to be the next knot location 
( u Ki v k)- At this point the algorithm proceeds to step b unless the maximum number 
of knot locations to be computed has been reached. 

As previously, we have used a QRP' decomposition of the coefficient matrix to solve the 
least squares problem for given knots and r value. This provides a stable and efficient means 
for solution of the problem with an indication if a matrix of less than full rank is encountered. 

In order to test the algorithms we have used a number of data sets. Several of these are 
based on previously published and widely available (x<y) datasets and parent functions. We 
have also used a few less readily available data sets that we are willing to share with anyone 
interested in obtaining them. Table 1 gives a summary of most of the data sets. 

As mentioned in the opening paragraph, the primary purpose of this paper is to discuss 
some further results we have obtained in allowing the coalescence of knots. A discussion of 
the algorithm we used is given in Section 2. Results of calculations obtained using several 
data sets and comparison with the previous method where only simple (but possibly near- 
repeated) knots are allowed, is given in Section 3. Section 4 summarizes our experiences. 
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n.m This refers to point set n and function m from [FR82], for n =1, 2, and 3, and 
m = 1 , * • • , 6 . n = 1 is 25 points, n = 2 is 33 points, n = 3 is 100 points, n = 4 
refers to the 200 point data set used in [MF92]. m=l is the humps and dip function, 
m = 2 is the cliff, m = 3 is the saddle, m = 4 is the gentle hill, m = 5 is the steep hill, 
and m = 6 is the sphere. In addition, m = 7 refers to the curved valley function from 
[NI78], 

GT This refers to the thinned glacier data consisting of 678 points, with certain contour 
lines removed, from [MF92]. 

GL This refers to the thinned glacier data consisting of 873 points. 

HF This is the data set from [MF92] generated to have point density approximately pro- 
portional to curvature, consisting of 500 points. 



Table 1: Data Sets Used Extensively in Tests 



2 Derivatives of Multiquadric Basis Functions and 
Multiple Knots 

Once again we have used Matlab 1 and LEASTSQ (a Levinburg-Marquardt method) to solve 
the nonlinear minimization problem 

N K 

min - ^2 a kB k {xi , y.) - c] 2 (4) 

k= 1 

where the minimization is taken over all (ufc,t^), r, the a^-, and c (with the last equation of 

(3) imposed as a constraint) as an equivalent two step minimization. Thus, for each given 
knot configuration and r value, the least squares solution of (3) is computed as a step toward 

(4) . This results in the solution of a simpler, but equivalent problem since 2 K parameters 
are eliminated from (4) by imposing the condition that the values of the aj and c be always 
taken as obtained from the least squares solution of (3). Hence, our final process is more 
properly written as 

N N 

min min ^[z,- - J2 a k B k (x t , y,) - c] 2 (5) 

t=l k= I 

where the inner minimization is over the aj and c (least squares solution of (3)), and the 
outer minimization is over the knot locations and the value of the parameter r. The global 
minimum of each of the two problems are clearly the same. Eq. (5) is the more restrictive, 
but any minimum of (4) is a local minimum of (5), else a better solution is attainable for (4). 

1 Math Works, 24 Prime Park Way, Natick. \1A 01760 
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This does not imply that the iterative methods employed to solve (5) would work equally 
well, nor find the same local minima, when applied to (4). 

When two knots approach each other, the system (3) becomes ill- conditioned, and the 
coefficients for the two knots tend to become large and of opposite sign. This reminds one of 
derivatives, in this case a directional derivative. Because the tw'o coefficients are not of the 
same magnitude, the appropriate replacement basis functions seem to be one of the basis 
functions, and a directional derivative of the basis function. In practice, the directional 
derivative is represented as a linear combination of the partial derivatives with respect to x 
and to y. To prove this is the proper change requires showing that the coefficients grow at no 
greater a rate than the reciprocal of the distance between the knots. The work of Narcowich 
and Ward [NW91] shows that the norm of the inverse of the multiquadric interpolation matrix 
grows no more rapidly than linear with the reciprocal of the minimum distance between data 
points. It seems likely a similar result holds for the norm of the pseudoinverse of the system 
matrix for (3). We have verified this computationally for some cases. 

The contribution of the two terms corresponding to the coalescent knots (say the jth 
and the kth) are of the form ajBj(x,y) + a^B^ix, y). The six parameters, a } , a^, u v 3 , 
Uk, and Vki in those terms are replaced by the five parameters in the terms a J B J (x,y) + 
bj(Bj(x,y)) x + Cj(Bj(x, 2/)) y . In the previous expression the subscripts x and y refer to 
partial differentiation. We will refer to such a knot as a double knot. 

Now, as it turns out, the coalescence of two knots, while the most prevalent, was by no 
means the only situation we encountered. The coalescence of three knots (or, once two knots 
have coalesced, the approach of a third knot to the location of a double knot) then results in 
the occurrence of a second derivative, a directional derivative of a directional derivative (the 
directions not necessarily being the same). In practice the replacement of the basis functions 
for the ith, jth, and kth knots is by the linear combination a 3 Bj + b J (B- l ) x + Cj(B-,) y + 
dj(Bj) xx -f ej(Bj) iy + fj(Bj) yy . While we have instances of near-coalescence of four knots, 
our present code does not attempt to handle such cases properly. 

The algorithm is required to change basis functions when it appears that knots are 
coalescing. We have used a simple idea that seems to be effective. We iterate the solution of 
(5) in stages having increasingly stringent convergence tolerances. As each tolerance level is 
met, a check for near- repeated knots is made, ’’near” being a changing tolerance consistent 
with the tolerance for the optimization. When pairs of knots are found within the tolerance, 
the new' basis is generated, and the optimization process continued. When no new repeated 
knots are obtained, the tolerance is decreased by a factor of two. We have observed that the 
inclusion of the derivative terms, vice the basis functions at near- repeated knots, generally 
yields essentially the same approximation, as it should, for both double and triple knots. The 
sequence of tolerance levels we used is: 0.02, 0.01, 0.005 for problems posed on the [0 , 1 ] 2 
square. 



3 Approximations with Repeated Knots 

The results of our investigation concern the ability of the algorithm to obtain good ap- 
proximations in comparison with the algorithm permitting only simple (but perhaps near- 
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repeated) knots, the degree of dependence on the initial guess for three different datasets 
(two in some detail), and a comparison of condition numbers of the coefficient matrices for 
the simple and multiple knot problems. 



dset 


r a 


7 'm 


rms s 


rms m 


grms a 


grms m 


knots. 


k TX ots fti 


3.1 


0.0976 


0.1588 


0.0081 


0.0059 


0.0115 


0.0107 


2nd, lnt 


Id, It, lnt 


3.2 


7e-6 


0.0002 


0.0074 


0.0073 


0.0186 


0.2477 


3nd,lnt 


3d,lt,lnq 


3.3 


0.3116 


0.2502 


0.0014 


0.0008 


0.0016 


0.0011 


lnt 


It 


3.4 


0.2885 


0.3118 
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0.0005 


0.0006 


0.0005 


2nd 


3d 


3.5 


0.3720 


0.3886 


0.0008 


0.0007 


0.0010 


0.0010 


3nd 


2d 


3.6 


2.3280 


2.3044 


0.0005 


0.0005 


0.0007 


0.0007 






3.7 


0.0990 


0.3304 


0.0403 


0.0279 


0.0733 


0.0324 


4nd 


4d 


4.1 


0.0958 


0.1588 


0.0070 


0.0017 


0.0071 


0.0024 


lnt 


It 


4.2 


2e-5 


2e-5 


0.0035 


0.0035 


0.0083 


0.0082 


4nd 


3d, lnd 


4.3 


0.2577 


0.2513 


0.0002 


0.0002 


0.0002 


0.0002 


4nd 


Id, lnd 


4.4 


0.4163 


0.4625 


7e-6 


6e-6 


7e-6 


6e-6 


lnd 


3d 


4.5 


0.3234 


0.3151 


0.0002 


0.0003 


0.0002 


0.0003 


4nd 


2nd 


4.6 


0.5668 


0.5680 


5e-5 


4e-5 


5e-5 


4e-5 


5nd 


2d,2nt,2nd 


4.7 


0.1111 


0.1134 


0.0115 


0.0114 


0.0202 


0.0287 


3nd,lnt 


5d,lt 



Table 2: Comparison of Results: Near- Repeated Knots vs. Multiple Knots 

We first describe the layout of the Table 2, which gives a comparison of the results of 
the simple knot and multiple knot algorithms for fourteen datasets. The results are for the 
surface fits obtained by starting from the initial guess obtained by the greedy algorithm 
[FH93], with r= 0.3, and closeness tolerance, ctol=0.1. The 3.n datasets were fit using twelve 
knot locations, while the 4.n datasets were fit using twenty knot locations. The columns 
of the table contain, in pairs, the values of r, the rms errors for the data points, the rms 
errors on a 21x21 grid, and information about the near-repeated or multiple knots obtained. 
The subscripts s and m denote the simple knot and multiple knot algorithms, respectively. 
The terms 2d and It, for example, refer to 2 double knots, and 1 triple knot. Likewise, lnd 
refers to 1 near-double knot. In some cases two double knots may have been obtained, with 
another simple knot nearby one of the double knots - such a case might be denoted by 2d, 
lnt. 

For a given dataset, the values of r obtained are mostly comparable, with larger values 
generally being obtained by the multiple knot algorithm. In the instances where the value 
of r is not larger, the values are essentially the same, the exception being dataset 3.3. The 
cases showing the most difficulty are the cases 3.2 and 4.2 (the cliff function), which has a 
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tendency toward small values of the parameter r. The effects of this tendency are discussed 
in a later paragraph. 

The rms errors are mostly comparable, with the errors obtained by the multiple knot 
algorithm generally smaller, with some ties, and one instance where the errors are bigger. 
This occurs for dataset 4.5, but the errors for both approximations are quite small for this 
case. 

The rms errors on the grid (grins) generally follow the same pattern, with one major 
exception. On the dataset 3.2, the errors for the grid are quite large. Closer inspection of 
this particular case reveals there is one grid point (near the triple, almost quadruple, knot) 
where the surface is badly ill-behaved. Coefficients of the first derivative terms are on the 
order of 5 to 20, while the coefficients of the two multiquadric basis functions are reasonably 
large, on the order of 10 3 . However, the value of r is very small in this case, about 1.5e-4, 
thus the derivatives of the basis functions very rapidly approach the value one in one of the 
coordinate directions (the value is about 0.99 at a distance of 0.001 from the knot). Hence, 
very rapid changes in the value of the function occur. While otherwise similar situations 
occur in other cases, the key difference here is the small value of r that occurs. 

In order to discover whether the initial guess was contributing to the problem, or whether 
the tendency toward a small r value and multiple knots was driven by the data, we ran a 
number initial guesses for the dataset 3.2. The initial conditions were obtained from the 
greedy algorithm with various values of r and closeness tolerance (ctol). The results of the 
computations are summarized in Table 3. The table gives the parameter ctol, the initial 
value of r, the final value of r, the rms errors at the data points and on the grid, and the 
number, multiplicity, and approximate location of the multiple knots. From this we see that 
while there is a tendency toward multiple knots, no particular locations are tended toward, 
although it seems natural that repeated knots invariably occur near the diagonal. The fact 
that the surface is a function of one variable (that is, y-x) with rapid changes occurring 
near the diagonal may account for this circumstance. The fact that the data is not regular 
may account for a large number of comparable local minima. One case is of particular note: 
The case in the second line involves three triple and one double knot, however two of the 
triple knots are close together, yielding a near-repeated knot of multiplicity six. This could 
be expected to lead to poor behavior of the surface. Closer inspection revealed that, even 
though the coefficients are rather large (order 10'), the surface is well behaved and yielded 
one of the better approximations among this group. This good behavior is undoubtedly 
enhanced by the fact that the r value is large for this surface. In the three cases where 
the r value is smaller than 10 -4 the resulting surface is poorly behaved near a multiple or 
near-multiple knot, this being reflected by the larger grms 0 values shown. 

We ran a number of different initial guesses for case 3.1 to determine both the robustness 
of the optimization routine and to confirm the tendency of the method to converge on 
multiple knots from a variety of initial guesses. Nine different initial guesses were made, 
six resulting from the greedy algorithm with different closeness tolerances and r values, and 
three with random initial knot locations. As is seen in Table 4, the results for the various 
initial guesses are strikingly similar: except for two cases the r values are all close to 0.16, all 
but one result in a triple or near-triple knot near the dip at about (.45, .78), and all but one 
have comparable accuracies for rms 0 and grms 0 . The one odd case had the starting value 
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of 0.01 for r. A wide-ranging difference is that different double and near-double knots may 
occur for the various runs. Unlike dataset 3.2, this case appears to be quite stable. 
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r 0 
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multiple knots 


0.1 


0.3 


0.0002 
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t@(.09,.14), 3d@(.38,.38, (.60, .60), (.89, .88)) 
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0.0239 


0.0060 


0.0085 
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0.0247 


0.0071 
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2d@(.38,.40), (.67, .66) 
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0.0113 


3d@(.14,.14),(.44,.44), (.68, .67) 
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0.0054 


0.0095 


4d@(.13,.ll),(.37,.39), (.69, .65), (.92, .91) 
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0.01 


5e-5 


0.0071 


0.0259 


d@(.69,.69), 2nd@(.14,.12), (.33, .40) 


0.2 


0.01 


0.0001 


0.0071 


0.0114 


d@(.44,.44) 


0.0 


0.01 


5e-5 


0.0071 


0.0260 


d@( .69,-69), 2nd@(.14,.ll), (.33, .40) 



Table 3: Results with Different Initial Guesses for Case 3.2 
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0.1 


0.1 


g 


0.1450 


0.0066 


0.0111 


t@(.43,.77) 


0.1 


0.01 


g 


0.0737 


0.0118 


0.0157 


d@(.72,.36) 


- 


0.3 


r 


0.1454 


0.0066 


0.0112 


t@(.43,.77) 


- 


0.3 


r 


0.1566 


0.0061 


0.0077 


t@(.42,.75) 


- 


0.3 


r 


0.1600 


0.0067 


0.0113 


t@(.43,.76), 2nd@(.83,.43), (.86, .23) 




0.3 


r 


0.1527 


0.0067 


0.0111 


nt@(.44,.75) 



Table 4: Results with Different Initial Guesses for Case 3.1 



We also ran some different initial conditions on dataset 3.7. The results in this case seem 
to indicate complications somewhere between the datasets 3.1 and 3.2. Dataset 3.7 has a 
tendency toward many multiple knots, perhaps because it has a more complex behavior of 
the surface than any of the other datasets. 
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The occurrence of near-multiple knots, as we noted previously, must have a unfavorable 
effect on the condition number of the coefficient matrix, and in order to investigate this we 
computed the condition number of the coefficient matrices for the fourteen cases of Table 2. 
A plot of the of the minimum separation distance divided by the r value vs. the condition 
number of the coefficient matrix is shown in Figure 1. There are 15 points for each of the 
simple knot and multiple knot algorithms. The first fourteen points correspond to the data 
in Table 2, in the order given, while the fifteenth is for the thinned glacier dataset, GL, 
using 25 knots. For the latter data the underlying function is unknown. The simple knot 
algorithm found 9 near-double knots. The multiple knot algorithm found four double knots 
and three triple knots. 

The points in Figure 1 for the simple knot algorithm have been labeled with their number. 
Because of the clustering of many of the points for the multiple knot algorithm, only those 
with condition number great than 10 5 have been labeled. 

What can be observed is that the condition numbers for the multiple knot algorithm are 
generally one to two orders of magnitude smaller than those for the simple knot algorithm. 
There are exceptions to this, most notably the datasets 3.1, 3.2, 3.7. and 4.5, where the be- 
havior is reversed, and 3.6 for which the r value is relatively large in both cases, contributing 
to the size of the condition number. 

Finally, for one of the parent surfaces, we include a plot of the underlying surface being 
approximated, along with the surfaces obtained from the simple knot algorithm and the 
multiple knot algorithm. The dataset is 3.7, the surface shown in Figure 2a being sampled 
at 100 points. Figures 2b and 2c, respectively, show the surface reconstructed using the 
simple knot and multiple knot algorithms with twelve knot points. In this case it is clear 
that the multiple knot algorithm has achieved a better fit. 



4 Conclusion and Ideas for Further Work 

We have formulated an algorithm to allow coalescence of near-repeated knots into multiple 
knots, bringing into the approximation the appropriate derivatives of the basis function in 
place of the basis functions for the nearby knots. As a general rule, this results in both 
better approximations, and more stable computations, although individual examples do not 
conform to this ideal. Our computations convince us that the occurrence of multiple knots 
is a natural phenomenon, and that algorithms for the problem need to account for the 
possibility. 

The tendency toward repeated knots leads to one possible idea for further exploration. 
This would be to begin initially with some (or all) knots being double knots to take advantage 
of the apparently greater fitting power of the multiquadric plus its directional derivative. We 
have not yet had the opportunity to pursue this idea. While the advantage of starting with 
multiple knots seems clear, the disadvantage is that coalescence of these knots skips odd- 
multiplicity knots, especially triple knots, which occur regularly in our computations. 

We have observed the tendency toward repeated knots with the LEASTSQ algorithm 
and another (FMINS - see [WOS5]) in the previous paper and feel the repeated knots are 
not occurring as an artifact of the optimization process. Nonetheless, an algorithm that took 
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advantage of the particular problem being solved could be more efficient and possibly avoid 
some of the local minima being found by the general minimization routines. 
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Figure Captions 



Figure 1: Minimum separation distance relative to r value vs. condition number for thirty 
problems. The o’s denote simple knot cases, numbered from 1-15, the first 14 being as 
ordered in Table 2. Number 15 is the thinned glacier data, GL, with 25 knots. The x’s 
denote the multiple knot cases - only those with condition number greater than 10 5 are 
labeled due to space considerations. 

Figure 2: a) The parent surface from which the data was obtained, b) The surface obtained 
with 12 knots using the simple knot algorithm, c) The surface obtained with 12 knots using 
the multiple knot algorithm. 
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